%matplotlib notebook
import numpy as np
import scipy.signal
import scipy.fft as sfft
import matplotlib.pylab as plt
from matplotlib import animation
from IPython.display import YouTubeVideo, HTML, Audio
from bokeh.layouts import column, row
from bokeh.models import CustomJS, ColumnDataSource, Slider
from bokeh.plotting import Figure, show, output_notebook
output_notebook()
2. Sistemas para el procesamiento de señales¶
2.1. Definición de sistema¶
Hasta ahora hemos realizado análisis de señales, es decir el estudio de las señales y sus propiedades en el dominio del tiempo y frecuencia
En esta unidad nos enfocaremos en el procesamiento de señales, es decir el diseño de sistemas que procesan una señal de entrada y producen una señal de salida
Usaremos
\(x[n]\) para denotar la señal (discreta) de entrada y \(X[k]\) su espectro
\(y[n]\) para denotar la señal (discreta) de salida e \(Y[k]\) su espectro
2.1.1. Ejemplos de sistemas¶
Utilizando sistemas podemos modificar una señal para mejorar su calidad o remover efectos indeseados
Un sistema para reducir el ruido de una señal de electroencefalograma (EEG)
Un sistema para mejorar una imagen fuera de foco (sharpening)
Un sistema para eliminar el eco de un audio
En lo que sigue realizaremos una clasificación de los sistemas según sus propiedades
2.2. Sistemas sin memoria¶
Diremos que un sistema \(\Phi\) es un sistema sin memoria si
es decir que la salida del sistema en un instante \(n\) dado depende solo de la entrada en ese mismo instante
Veamos algunos ejemplos
2.2.1. Sistema atenuador/amplificador ideal¶
donde \(A>0\) se llama ganancia
Este sistema puede atenuar la entrada si \(0<A<1\) y amplificar si \(A>1\)
2.2.2. Sistema saturador (clamp)¶
Este sistema limita los valores de la señal de entrada en un rango fijo
2.2.3. Sistema rectificador¶
Este sistema eliminar la parte negativa de la señal de entrada
2.3. Sistema Lineal¶
Diremos que un sistema \(\Phi\) es lineal si cumple con las siguientes propiedades
2.3.1. Homogeneidad¶
Un cambio en la amplitud de la entrada produce un cambio equivalente en la salida
2.3.2. Aditividad¶
Señales que se suman en la entrada producen señales que se suman en la salida
Es decir que las señales pasan por el sistema sin interactuar entre ellas
2.3.3. Otras propiedades de los sistemas lineales¶
Producto de las propiedades anteriores se tiene que una cascada de sistemas lineales forma un sistema lineal equivalente
Y la cascada de sistemas es conmutativa: El orden de los sistemas en la cascada no altera el resultado final
Los sistemas lineales también cumplen el Principio de superposición
Si descomponemos una señal en \(M\) componentes: \(x[n] = x_1[n] + x_2[n] + \ldots + x_M[n]\)
Y aplicamos un sistema lineal a cada componente: \(y_j[n] = \Phi(x_j[n])\)
Podemos recuperar la salida total usando aditividad: \(y_1[n] + y_2[n] + \ldots + y_M[n] = y[n]\)
2.4. Sistemas con memoria¶
Un sistema \(\Phi\) es un sistema con memoria si su salida actual depende sólo de la entrada actual, las entradas anteriores o las salidas anteriores
esto también se conoce como sistema causal
Un sistema con memoria no-causal usa entradas futuras (es decir \(x[n+1]\), \(x[n+2]\), …) y por ende solo se puede implementar de forma offline, es decir una vez que sea ha observado toda la señal
Veamos algunos ejemplos de sistemas con memoria
2.4.1. Sistema con un retardo (delay)¶
Definido como
donde
la salida depende solo de «una» entrada anterior
el valor de m define que tan «antigua» es la entrada pasada
El delay o retarno no afecta la amplitud de los componentes frecuenciales de la señal pero si su fase, como muestra la siguiente animación
%%capture
fig, ax = plt.subplots(3, figsize=(6, 6), tight_layout=True)
n = np.arange(0, 200, step=1)
x = lambda m: np.sin(2.0*np.pi*0.05*(n-m))
f = sfft.fftshift(sfft.fftfreq(d=1, n=len(n)))
def update(m):
ax[0].cla(); ax[0].plot(n, x(m));
X = sfft.fftshift(sfft.fft(x(m)))
Xm = np.absolute(X); Xp = np.angle(X)
# Espectro de magnitud:
ax[1].cla(); ax[1].plot(f, Xm);
# Espectro de fase enmascarado con el espectro de magnitud
ax[2].cla(); ax[2].plot(f, Xm*Xp/np.amax(Xm));
ax[2].set_ylim([-np.pi, np.pi])
angle_delay = Xp[np.argmax(Xm)]
ax[2].set_title("%0.4f [rad], %0.0f [deg]" % (angle_delay, 180*angle_delay/np.pi))
return ()
anim = animation.FuncAnimation(fig, update, frames=40, interval=150, blit=True)
HTML(anim.to_html5_video())
2.4.2. Sistema reverberador o eco¶
Definido como
donde
la salida depende de una entrada «pasada» y la entrada actual
la ganancia controla si el «eco» se atenua o amplifica
Al contrario del sistema anterior, el eco si puede modificar el espectro de amplitud.
Notemos el efecto de interferencia constructiva y destructiva al modificar el retardo, como muestra la siguiente animación
%%capture
fig, ax = plt.subplots(3, figsize=(6, 6), tight_layout=True)
n = np.arange(0, 200, step=1)
x = lambda m: np.sin(2.0*np.pi*0.05*(n-m))
f = sfft.fftshift(sfft.fftfreq(d=1, n=len(n)))
A = 1.
def update(m):
y = x(0) + A*x(m)
X = sfft.fftshift(sfft.fft(y))
for ax_ in ax:
ax_.cla()
ax[0].plot(n, x(0), label='x[n]')
ax[0].plot(n, A*x(m), label='A*x[n-m]')
ax[0].legend()
ax[1].plot(n, y);
ax[1].set_ylim([-A-1.1, A+1.1])
ax[2].plot(f, np.absolute(X));
ax[2].set_ylim([-20, (1+A)*len(n)/2 + 20])
anim = animation.FuncAnimation(fig, update, frames=40, interval=150, blit=False)
HTML(anim.to_html5_video())
Ejemplo de interferencia destructiva en una onda mecánica: https://www.youtube.com/watch?v=IU8xeJlJ0mk
2.4.3. Sistemas con múltiples ecos¶
Pueden combinarse más retardos para hacer un sistema reverberante más complejo
Por ejemplo
Que se escucha como
Fs = 44100;
n = np.arange(0, 4, step=1.0/Fs)
x = lambda m: np.sin(2.0*np.pi*880*(n-m))*np.exp(-(n-m)**2/0.5**2)*np.heaviside(n-m, 0)
y = x(0) + 0.5*x(1.) + 0.25*x(2.) + 0.125*x(3.)
p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(n, y, line_width=3)
p1.xaxis[0].axis_label = 'Tiempo [s]'
show(p1)
Audio(y, rate=Fs, normalize=False)
2.5. Sistema FIR¶
Generalizando el ejemplo de sistema lineal reverberante a \(L\) retardos llegamos a
que se puede modelar como una convolución discreta entre \(h\) y \(x\)
Este sistema se conoce como
Sistema FIR (finite impulse response)
Sistema MA (moving average)
Sistema todo-zeros
y es de orden L (posee L+1 coeficientes)
2.5.1. Intepretación como media movil (MA)¶
El sistema FIR es equivalente a una media movil ponderada que se aplica sobre la entrada
Los coeficientes del sistema son los ponderadores
Por ejemplo sea un sistema de 3 coeficientes \(h[0]=a\), \(h[1]=b\) y \(h[2]=c\)
donde cada salida se calcula a partir de
Para obtener el valor de \(y[0]\) e \(y[1]\) se deben establecer «condiciones de borde», como por ejemplo \(x[-2] = x[-1]= 0\)
2.5.2. Ejemplo: Eliminando ruido blanco aditivo¶
Si tenemos una señal contaminada con ruido blanco aditivo podemos usar un sistema FIR promediador para «suavizar la contaminación»
La animación muestra un sistema FIR con \(L\) coeficientes idénticos e iguales a \(1/L\) que se convoluciona con la señal contaminada intentando recuperar la «señal limpia»
%%capture
np.random.seed(0)
n = np.arange(0, 100, step=1)
C = 5*np.exp(-0.5*(n[:, np.newaxis] - n[:, np.newaxis].T)**2/10**2)
# Señal de entrada original
x_clean = np.random.multivariate_normal(np.zeros_like(n), C)
# Señal de entrada contaminada
x = x_clean + 2*np.random.randn(len(n))
# Sistema
L = 10; h = np.ones(shape=(L,))/L;
# Acumulator
y = np.zeros_like(n, dtype=np.float)
fig, ax = plt.subplots(3, figsize=(7, 6), tight_layout=True)
ax[0].plot(n, x, 'k.');
ylims = ax[0].get_ylim()
def init():
global y
y = np.zeros_like(n, dtype=np.float)
return ()
def update(m):
ax[1].cla(); ax[2].cla()
c = np.zeros_like(n, dtype=np.float);
c[m:m+L] = h
ax[1].plot(n, c);
y[m] = np.sum(h*x[m:m+L])
ax[2].plot(n, y, label='señal recuperada');
ax[2].plot(n, x_clean, 'g', label='señal limpia');
ax[2].legend()
ax[2].set_ylim(ylims)
ax[2].plot([m, m], [ylims[0], ylims[1]], 'r--', alpha=0.5)
return ()
anim = animation.FuncAnimation(fig, update,init_func=init,
frames=100-len(h), interval=150, blit=True)
HTML(anim.to_html5_video())
2.5.3. Ejemplo: Encontrando cambios en una señal¶
Podemos usar un sistema «diferenciador» para detectar cuando una señal cambia de valor como se muestra en la siguiente animación
%%capture
fig, ax = plt.subplots(3, figsize=(7, 6), tight_layout=True)
n = np.arange(0, 100, step=1)
x = np.zeros_like(n, dtype=np.float)
x[20:] += 1.; x[40:] += 1.; x[80:] -= 1.;
ax[0].plot(n, x)
# System:
h = np.array([-0.5, 0.5])
# Acumulator
y = np.zeros_like(n, dtype=np.float)
def init():
global y
y = np.zeros_like(n, dtype=np.float)
return ()
def update(m):
c = np.zeros_like(n, dtype=np.float); c[m:m+len(h)] = h
ax[1].cla(); ax[1].plot(n, c);
y[m] = np.sum(h*x[m:m+len(h)])
ax[2].cla(); ax[2].plot(n, y);
ax[2].plot([m, m], [-0.5, 0.5], 'r--', alpha=0.5)
return ()
anim = animation.FuncAnimation(fig, update, init_func=init,
frames=100-len(h), interval=150, blit=True)
HTML(anim.to_html5_video())
2.5.4. Ejemplo: Removiendo una tendencia¶
En el ejemplo siguiente tenemos una señal sinusoidal que nos interesa separar de una señal más lenta/suave denominada tendencia
Podemos usar un sistema FIR para eliminar la tendencia
%%capture
fig, ax = plt.subplots(3, figsize=(7, 6), tight_layout=True)
np.random.seed(0);
n = np.arange(0, 100, step=1)
C = np.exp(-0.5*(n[:, np.newaxis] - n[:, np.newaxis].T)**2/30**2)
x = np.sin(2.0*np.pi*0.1*n) + 5*np.random.multivariate_normal(np.zeros_like(n), C)
ax[0].plot(n, x); ylims = ax[0].get_ylim()
# System:
L = 5; h = -np.ones(shape=(L,))/L; h[L//2] += 1
# Acumulator
y = np.zeros_like(n, dtype=np.float)
def init():
global y
y = np.zeros_like(n, dtype=np.float)
return ()
def update(m):
c = np.zeros_like(n, dtype=np.float); c[m:m+len(h)] = h
ax[1].cla(); ax[1].plot(n, c);
y[m] = np.sum(h*x[m:m+len(h)])
ax[2].cla(); ax[2].plot(n, y);
ax[2].plot([m, m], [-0.5, 0.5], 'r--', alpha=0.5)
return ()
anim = animation.FuncAnimation(fig, update, init_func=init,
frames=100-len(h), interval=200, blit=True)
HTML(anim.to_html5_video())
2.6. Convolución con scipy¶
En la práctica, si queremos convolucionar una señal con nuestro sistema, podemos usar la función de scipy
scipy.signal.convolve(in1, # Señal de entrada
in2, # Coeficientes del sistema
mode='full',
method='auto'
)
donde el argumento method puede ser
direct: Realiza la convolución en el dominio del tiempofft: Realiza la convolución multiplicando los espectrosauto: Se decide automaticamente en base al largo de las señales
y el argumento mode indica donde se hace la convolución
Sea una señal \(x=[a,b,c]\) y un sistema \(h=[d, e]\)
Si uso
mode=validel resultado será \(y=[ad+be, bd+ce]\)Si uso
mode=sameel resultado será \(y=[ae, ad+be, bd+ce]\), es decir se agregan ceros al principio de \(x\) tal que \(y\) sea del mismo largo que \(x\)Si uso
mode=fullel resultado será \(y=[ae, ad+be, bd+ce, cd]\), es decir se agregan ceros al principio y al final de \(x\)
2.6.1. Ejemplo: Eliminando ruido v 2.0¶
Esta vez usamos scipy.signal.convolve y probamos dos sistemas FIR
coeficientes idénticos e iguales a \(1/L\) (rectangular)
coeficientes que decaen suavemente a cero (hamming)
from bokeh.palettes import Dark2_5 as palette
np.random.seed(0)
n = np.arange(0, 100, step=1)
C = 5*np.exp(-0.5*(n[:, np.newaxis] - n[:, np.newaxis].T)**2/10**2)
# Señal de entrada original
x_clean = np.random.multivariate_normal(np.zeros_like(n), C)
# Señal de entrada contaminada
x = x_clean + 2*np.random.randn(len(n))
L = 20;
p = [Figure(plot_width=600, plot_height=230, toolbar_location="below") for k in range(3)]
p[0].scatter(n, x, color='black')
p[2].line(n, x_clean, line_width=2, legend_label = 'señal limpia');
for name, h, color in zip(['rect', 'hamming'],
[np.ones(shape=(L,)), scipy.signal.hamming(L)],
palette):
h = h/np.sum(h)
p[1].line(np.arange(L), h, color=color, legend_label=name, line_width=2)
p[2].line(n, scipy.signal.convolve(x, h, mode='same', method='auto'),
color=color, legend_label=name, line_width=2);
p[2].legend.location = "top_left"
show(column(p))